Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
Produced with R version 3.2.3 and pomp version 1.3.1.2.
To discuss forecasting and fitted values, both as ways to assess model fit and as potentially the major goal of time series analysis.
To teach some forecasting methods based on POMP models.
To demonstrate a forecasting analysis using POMP methods and the pomp computational environment.
Given data up to time \(t_n\), and a fitted model, what is our prediction for the future values of a process \(X_{{n+1}}, X_{n+1},\dots\)?
For a POMP model, the answer is fairly clear. According to the structure of a POMP model, the data up to time \(t_n\) affect our estimate of the latent variable \(X_n\), but do not otherwise affect the future evolution of the system.
If we have constructed the filtering distribution, \[ f_{X_n|Y_{1:n}}(x_n|{y_{1:n}^*}{\, ; \,}\hat\theta),\] for some appropriate parameter estimat \(\hat\theta\), our forecast distribution for \(X_{n+p}\) given \({y_{1:n}^*}\) involves combining the filtering distribution with the transition density for \(X_{n+p}\) given \(X_n\). We can write this as
\[\begin{eqnarray} f_{X_{n+p}|Y_{1:n}}(x_{n+p}|{y_{1:n}^*}{\, ; \,}\hat\theta) &=& \int f_{X_n|Y_{1:n}}(x_n|{y_{1:n}^*}{\, ; \,}\hat\theta) f_{X_{n+p}|X_n}(x_{n+p}{{\, | \,}}x_n{\, ; \,}\hat\theta) \, dx_{n+p} \\ &=& \int f_{X_n|Y_{1:n}}(x_n|{y_{1:n}^*}{\, ; \,}\hat\theta) \prod_{k=1}^p f_{X_{n+k}|X_{n+k-1}}(x_{n+k}{{\, | \,}}x_{n+k-1}{\, ; \,}\hat\theta) \, dx_{n+1}\dots dx_{n+p} \end{eqnarray}\]
A closely related task is to forecast \(Y_{n+p}\) rather than \(X_{n+p}\).
Since an ARIMA model is a POMP model, this theory also applies to ARIMA analysis and can be implemented in R by applying predict to the output of an arima.
For a nonlinear POMP analysis in pomp, the Monte Carlo approach is to represent the filtering distribution by \(J\) particles produced from pfilter. Then, the Monte Carlo representation of \(f_{X_{n+p}|Y_{1:n}}(x_{n+p}|{y_{1:n}^*}{\, ; \,}\hat\theta)\) comes from applying rprocess to these filter particles.
We may wish to extend this approach to generate a forecast that takes into account uncertainty in \(\theta\). Weighting a set of candidate parameter values by their likelihood given the data provides one way to do thatโan empirical Bayes approach.
The scenario. Letโs situate ourselves at the beginning of October 2014. Ebola has quite suddenly turned into a threat to global health and global trade. The WHO situation report contained data on the number of Ebola cases in each of Guinea, Sierra Leone, and Liberia. Key questions included:
Download the data from the WHO Situation Report of 1 October 2014:
base_url <- "http://kingaa.github.io/sbied/"
read.csv(paste0(base_url,"data/ebola_data.csv"),stringsAsFactors=FALSE,
colClasses=c(date="Date")) -> dat
sapply(dat,class)
## week date country cases
## "integer" "Date" "character" "numeric"
head(dat)
## week date country cases
## 1 1 2014-01-05 Guinea 2.244
## 2 2 2014-01-12 Guinea 2.244
## 3 3 2014-01-19 Guinea 0.073
## 4 4 2014-01-26 Guinea 5.717
## 5 5 2014-02-02 Guinea 3.954
## 6 6 2014-02-09 Guinea 5.444
Supplementing these data are population estimates for the three countries.
## Population sizes in Guinea, Liberia, and Sierra Leone (census 2014)
populations <- c(Guinea=10628972,Liberia=4092310,SierraLeone=6190280)
dat %>%
ggplot(aes(x=date,y=cases,group=country,color=country))+
geom_line()
Many of the early modeling efforts used variants on the simple SEIR model. Here, weโll focus on a variant that attempts a more accurate description of the duration of the latent period. Specifically, this model assumes that the amount of time an infection remains latent is \[\mathrm{LP} \sim \mathrm{Gamma}\left(m,\frac{1}{m\,\alpha}\right),\] where \(m\) is an integer. This means that the latent period has expectation \(1/\alpha\) and variance \(1/(m\,\alpha)\). In this document, weโll fix \(m=3\).
We implement Gamma distributions using the so-called linear chain trick.
rSim <- Csnippet('
double lambda, beta;
double *E = &E1;
beta = R0 * gamma; // Transmission rate
lambda = beta * I / N; // Force of infection
int i;
// Transitions
// From class S
double transS = rbinom(S, 1.0 - exp(- lambda * dt)); // No of infections
// From class E
double transE[nstageE]; // No of transitions between classes E
for(i = 0; i < nstageE; i++){
transE[i] = rbinom(E[i], 1.0 - exp(-nstageE * alpha * dt));
}
// From class I
double transI = rbinom(I, 1.0 - exp(-gamma * dt)); // No of transitions I->R
// Balance the equations
S -= transS;
E[0] += transS - transE[0];
for(i=1; i < nstageE; i++) {
E[i] += transE[i-1] - transE[i];
}
I += transE[nstageE-1] - transI;
R += transI;
N_EI += transE[nstageE-1]; // No of transitions from E to I
N_IR += transI; // No of transitions from I to R
')
The deterministic skeleton is an ODE.
skel <- Csnippet('
double lambda, beta;
const double *E = &E1;
double *DE = &DE1;
beta = R0 * gamma; // Transmission rate
lambda = beta * I / N; // Force of infection
int i;
// Balance the equations
DS = - lambda * S;
DE[0] = lambda * S - nstageE * alpha * E[0];
for (i=1; i < nstageE; i++)
DE[i] = nstageE * alpha * (E[i-1]-E[i]);
DI = nstageE * alpha * E[nstageE-1] - gamma * I;
DR = gamma * I;
DN_EI = nstageE * alpha * E[nstageE-1];
DN_IR = gamma * I;
')
\(C_t | H_t\) is negative binomial with \(\expect{C_t|H_t} = \rho\,H_t\) and \({\mathrm{Var}}{C_t|H_t} = \rho\,H_t\,(1+k\,\rho\,H_t)\).
dObs <- Csnippet('
double f;
if (k > 0.0)
f = dnbinom_mu(nearbyint(cases),1.0/k,rho*N_EI,1);
else
f = dpois(nearbyint(cases),rho*N_EI,1);
lik = (give_log) ? f : exp(f);
')
rObs <- Csnippet('
if (k > 0) {
cases = rnbinom_mu(1.0/k,rho*N_EI);
} else {
cases = rpois(rho*N_EI);
}')
toEst <- Csnippet('
const double *IC = &S_0;
double *TIC = &TS_0;
TR0 = log(R0);
Trho = logit(rho);
Tk = log(k);
to_log_barycentric(TIC,IC,4);
')
fromEst <- Csnippet('
const double *IC = &S_0;
double *TIC = &TS_0;
TR0 = exp(R0);
Trho = expit(rho);
Tk = exp(k);
from_log_barycentric(TIC,IC,4);
')
The following function constructs a pomp object to hold the data for any one of the countries.
ebolaModel <- function (country=c("Guinea", "SierraLeone", "Liberia"),
timestep = 0.1, nstageE = 3) {
ctry <- match.arg(country)
pop <- unname(populations[ctry])
nstageE <- as.integer(nstageE)
globs <- paste0("static int nstageE = ",nstageE,";")
dat <- subset(dat,country==ctry,select=-country)
## Create the pomp object
dat %>%
extract(c("week","cases")) %>%
pomp(
times="week",
t0=min(dat$week)-1,
globals=globs,
statenames=c("S","E1","I","R","N_EI","N_IR"),
zeronames=c("N_EI","N_IR"),
paramnames=c("N","R0","alpha","gamma","rho","k",
"S_0","E_0","I_0","R_0"),
nstageE=nstageE,
dmeasure=dObs, rmeasure=rObs,
rprocess=discrete.time.sim(step.fun=rSim, delta.t=timestep),
skeleton=skel, skeleton.type="vectorfield",
toEstimationScale=toEst,
fromEstimationScale=fromEst,
initializer=function (params, t0, nstageE, ...) {
all.state.names <- c("S",paste0("E",1:nstageE),"I","R","N_EI","N_IR")
comp.names <- c("S",paste0("E",1:nstageE),"I","R")
x0 <- setNames(numeric(length(all.state.names)),all.state.names)
frac <- c(params["S_0"],rep(params["E_0"]/nstageE,nstageE),params["I_0"],params["R_0"])
x0[comp.names] <- round(params["N"]*frac/sum(frac))
x0
}
) -> po
}
ebolaModel("Guinea") -> gin
ebolaModel("SierraLeone") -> sle
ebolaModel("Liberia") -> lbr
King et al. (2015) estimated parameters for this model for each country. A large Latin hypercube design was used to initiate a large number of iterated filtering runs. Profile likelihoods were computed for each country against the parameters \(k\) (the measurement model overdispersion) and \(R_0\) (the basic reproductive ratio). Full details are given on the datadryad.org site. The following loads the results of these calculations.
options(stringsAsFactors=FALSE)
profs <- read.csv(paste0(base_url,"/ebola/ebola-profiles.csv"))
The following plots the profile likelihoods. The horizontal line represents the critical value of the likelihood ratio test for \(p=0.01\).
require(reshape2)
require(plyr)
require(magrittr)
require(ggplot2)
theme_set(theme_bw())
profs %>%
melt(id=c("profile","country","loglik")) %>%
subset(variable==profile) %>%
ddply(~country,mutate,dll=loglik-max(loglik)) %>%
ddply(~country+profile+value,subset,loglik==max(loglik)) %>%
ggplot(mapping=aes(x=value,y=dll))+
geom_point(color='red')+
geom_hline(yintercept=-0.5*qchisq(p=0.99,df=1))+
facet_grid(country~profile,scales='free')+
labs(y=expression(l))
Parameter estimation is the process of finding the parameters that are โbestโ, in some sense, for a given model, from among the set of those that make sense for that model. Model selection, likewise, aims at identifying the โbestโ model, in some sense, from among a set of candidates. One can do both of these things more or less well, but no matter how carefully they are done, the best of a bad set of models is still bad.
Letโs investigate the model here, at its maximum-likelihood parameters, to see if we can identify problems. The guiding principle in this is that, if the model is โgoodโ, then the data are a plausible realization of that model. Therefore, we can compare the data directly against model simulations. Moreover, we can quantify the agreement between simulations and data in any way we like. Any statistic, or set of statistics, that can be applied to the data can also be applied to simulations. Shortcomings of the model should manifest themselves as discrepancies between the model-predicted distribution of such statistics and their value on the data.
pomp provides tools to facilitate this process. Specifically, the probe function applies a set of user-specified probes or summary statistics, to the model and the data, and quantifies the degree of disagreement in several ways.
Letโs see how this is done using the model for the Guinean outbreak.
library(pomp)
library(plyr)
library(reshape2)
library(magrittr)
options(stringsAsFactors=FALSE)
profs %>%
subset(country=="Guinea") %>%
subset(loglik==max(loglik),
select=-c(loglik,loglik.se,country,profile)) %>%
unlist() -> coef(gin)
simulate(gin,nsim=20,as.data.frame=TRUE,include.data=TRUE) %>%
mutate(date=min(dat$date)+7*(time-1),
is.data=ifelse(sim=="data","yes","no")) %>%
ggplot(aes(x=date,y=cases,group=sim,color=is.data,
alpha=is.data))+
geom_line()+
guides(color=FALSE,alpha=FALSE)+
scale_color_manual(values=c(no=gray(0.6),yes='red'))+
scale_alpha_manual(values=c(no=0.5,yes=1))
The simulations appear to be growing a bit more quickly than the data. Letโs try to quantify this. First, weโll write a function that estimates the exponential growth rate by linear regression. Then, weโll apply it to the data and to 500 simulations.
growth.rate <- function (y) {
cases <- y["cases",]
fit <- lm(log1p(cases)~seq_along(cases))
unname(coef(fit)[2])
}
probe(gin,probes=list(r=growth.rate),nsim=500) %>% plot()
Do these results bear out our suspicion that the model and data differ in terms of growth rate?
The simulations also appear to be more highly variable around the trend than do the data.
growth.rate.plus <- function (y) {
cases <- y["cases",]
fit <- lm(log1p(cases)~seq_along(cases))
c(r=unname(coef(fit)[2]),sd=sd(residuals(fit)))
}
probe(gin,probes=list(growth.rate.plus),
nsim=500) %>% plot()
Letโs also look more carefully at the distribution of values about the trend using the 1st and 3rd quantiles. Also, it looks like the data are less jagged than the simulations. We can quantify this using the autocorrelation function (ACF).
log1p.detrend <- function (y) {
cases <- y["cases",]
y["cases",] <- as.numeric(residuals(lm(log1p(cases)~seq_along(cases))))
y
}
probe(gin,probes=list(
growth.rate.plus,
probe.quantile(var="cases",prob=c(0.25,0.75)),
probe.acf(var="cases",lags=c(1,2,3),type="correlation",
transform=log1p.detrend)
),nsim=500) %>% plot()
Apply probes to investigate the extent to which the model is an adequate description of the data from the Sierra Leone outbreak. Have a look at the probes provided with pomp: ?basic.probes. Try also to come up with some informative probes of your own. Discuss the implications of your findings.
Up to now, weโve primarily focused on using POMP models to answer scientific questions. Of course, we can also use them to make forecasts. The key issues are to do with quantifying the forecast uncertainty. This arises from four sources:
Here, weโll explore how we can account for the first three of these in making forecasts for the Sierra Leone outbreak.
require(pomp)
require(plyr)
require(reshape2)
require(magrittr)
options(stringsAsFactors=FALSE)
set.seed(988077383L)
## forecast horizon
horizon <- 13
profs %>%
subset(country=="SierraLeone") %>%
subset(loglik==max(loglik),
select=-c(loglik,loglik.se,country,profile)) %>%
unlist() -> mle
## Weighted quantile function
wquant <- function (x, weights, probs = c(0.025,0.5,0.975)) {
idx <- order(x)
x <- x[idx]
weights <- weights[idx]
w <- cumsum(weights)/sum(weights)
rval <- approx(w,x,probs,rule=1)
rval$y
}
profs %>%
subset(country=="SierraLeone",
select=-c(country,profile,loglik.se)) %>%
subset(loglik>max(loglik)-0.5*qchisq(df=1,p=0.99)) %>%
melt(variable.name="parameter") %>%
ddply(~parameter,summarize,
min=min(value),max=max(value)) %>%
subset(parameter!="loglik") %>%
melt(measure=c("min","max")) %>%
acast(parameter~variable) -> ranges
params <- sobolDesign(lower=ranges[,'min'],
upper=ranges[,'max'],
nseq=20)
plot(params)
require(foreach)
require(doMC)
require(iterators)
registerDoMC(cores=4)
set.seed(887851050L,kind="L'Ecuyer")
foreach(p=iter(params,by='row'),
.inorder=FALSE,
.combine=rbind,
.options.multicore=list(preschedule=TRUE,set.seed=TRUE)
) %dopar%
{
require(pomp)
M1 <- ebolaModel("SierraLeone")
pf <- pfilter(M1,params=unlist(p),Np=2000,save.states=TRUE)
pf$saved.states %>% tail(1) %>% melt() %>%
dcast(rep~variable,value.var="value") %>%
ddply(~rep,summarize,S_0=S,E_0=E1+E2+E3,I_0=I,R_0=R) %>%
melt(id="rep") %>% acast(variable~rep) -> x
pp <- parmat(unlist(p),ncol(x))
simulate(M1,params=pp,obs=TRUE) %>%
melt() %>%
mutate(time=time(M1)[time],
period="calibration",
loglik=logLik(pf)) -> calib
M2 <- M1
time(M2) <- max(time(M1))+seq_len(horizon)
timezero(M2) <- max(time(M1))
pp[rownames(x),] <- x
simulate(M2,params=pp,obs=TRUE) %>%
melt() %>%
mutate(time=time(M2)[time],
period="projection",
loglik=logLik(pf)) -> proj
rbind(calib,proj)
} %>% subset(variable=="cases",select=-variable) %>%
mutate(weight=exp(loglik-mean(loglik))) %>%
arrange(time,rep) -> sims
ess <- with(subset(sims,time==max(time)),weight/sum(weight))
ess <- 1/sum(ess^2); ess
## [1] 10126.79
sims %>% ddply(~time+period,summarize,prob=c(0.025,0.5,0.975),
quantile=wquant(value,weights=weight,probs=prob)) %>%
mutate(prob=mapvalues(prob,from=c(0.025,0.5,0.975),
to=c("lower","median","upper"))) %>%
dcast(period+time~prob,value.var='quantile') %>%
mutate(date=min(dat$date)+7*(time-1)) -> simq
simq %>% ggplot(aes(x=date))+
geom_ribbon(aes(ymin=lower,ymax=upper,fill=period),alpha=0.3,color=NA)+
geom_line(aes(y=median,color=period))+
geom_point(data=subset(dat,country=="SierraLeone"),
mapping=aes(x=date,y=cases),color='black')+
labs(y="cases")
King, A. A., M. Domenech de Cellรจs, F. M. G. Magpantay, and P. Rohani. 2015. Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to ebola. Proc. R. Soc. Lond. B 282:20150347.